BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, Float64[], tspan, Float64[]) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true)
Parsing parameters...done
Creating parameters...done
Parsing species...done
Creating species...done
Creating species and parameters for evaluating expressions...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 414s / 98.9% 69.8GiB / 99.5%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
ODEProb DenseJac 1 289s 70.5% 289s 45.9GiB 66.1% 45.9G
iB
ODEProb SparseJac 1 61.9s 15.1% 61.9s 14.7GiB 21.2% 14.7G
iB
ODEProb No Jac 1 33.2s 8.11% 33.2s 4.76GiB 6.85% 4.76G
iB
Create ODESys 1 14.0s 3.43% 14.0s 2.62GiB 3.77% 2.62G
iB
Parse Network 1 11.5s 2.81% 11.5s 1.42GiB 2.05% 1.42G
iB
──────────────────────────────────────────────────────────────────────────
──
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show length(parameters(rn)) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 length(parameters(rn)) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = ModelingToolkit.varmap_to_vars(nothing, species(rn); defaults=ModelingToolkit.defaults(rn)) du = copy(u) p = ModelingToolkit.varmap_to_vars(nothing, parameters(rn); defaults=ModelingToolkit.defaults(rn)) @timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.) @timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.) densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.)
We also time the ODE rhs function with BenchmarkTools as it is more accurate given how fast evaluating f is:
@btime oprob.f($du,$u,$p,0.)
37.820 μs (3 allocations: 512 bytes)
Now we time the Jacobian functions, including compilation time in the first evaluations
J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 807s / 97.4% 86.6GiB / 99.3%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
SparseJac Eval1 1 306s 38.9% 306s 8.68GiB 10.1% 8.68G
iB
ODEProb DenseJac 1 289s 36.8% 289s 45.9GiB 53.4% 45.9G
iB
ODE rhs Eval1 1 68.0s 8.66% 68.0s 5.88GiB 6.84% 5.88G
iB
ODEProb SparseJac 1 61.9s 7.89% 61.9s 14.7GiB 17.1% 14.7G
iB
ODEProb No Jac 1 33.2s 4.23% 33.2s 4.76GiB 5.53% 4.76G
iB
Create ODESys 1 14.0s 1.79% 14.0s 2.62GiB 3.04% 2.62G
iB
Parse Network 1 11.5s 1.47% 11.5s 1.42GiB 1.65% 1.42G
iB
DenseJac Eval1 1 2.28s 0.29% 2.28s 2.01GiB 2.34% 2.01G
iB
SparseJac Eval2 1 92.3μs 0.00% 92.3μs 944B 0.00% 94
4B
ODE rhs Eval2 1 76.3μs 0.00% 76.3μs 688B 0.00% 68
8B
──────────────────────────────────────────────────────────────────────────
──
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)
